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difference  method.  The  resultin'  set  of  algebraic  equations  were 
solved  usin'  the  C-auss-Seidel  iterative  technique.  Tne  temperatures 
at  the  nodal  ;x>ints  we re  substituted  into  a  numerical  approximation  of 
the  variational  integral.  The  variational  integral  approximation  was 
used  to  determine  when  to  stop  the  iterative  process.  Tne  variational 
integral  stopping  criterion  was  compared  to  a  stoopin'  criterion  that 
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I  Introduction 


Background 

In  practice,  the  number  of  boundary  value  problems  that  can  be 
solved  analytically  is  quite  limited.  Consequently,  numerical  techni¬ 
ques  are  used  to  approximate  the  partial  differential  equations  associ¬ 
ated  with  the  boundary  value  problems.  These  numerical  techniques 
often  lead  to  large  sets  of  simultaneous  equations.  These  ec  _ions 
are  usually  sparse  (they  contain  few  non-zero  elements  in  the  -coeffi¬ 
cient  matrix) . 

Iterative  techniques  are  often  used  to  solve  large  sets  of  simul¬ 
taneous  equations.  These  techniques  have  two  major  advantages  over 
methods  that  solve  the  simultaneous  equations  directly.  First,  itera¬ 
tive  methods  avoid  calculations  with  the  many  zero  elements  in  the 
coefficient  matrix  while  direct  methods  perform  calculations  on  every 
element  of  the  coefficient  matrix.  Secondly,  iterative  methods  are  not 
as  susceptible  to  the  accumulation  of  round-off  error. 

An  important  difficulty  with  the  iterative  techniques  is 
determining  when  a  sufficient  number  of  iterations  has  been  completed 
to  ensure  an  accurate  solution.  This  thesis  addresses  the  problem  of 
when  to  step  iterating. 

Many  boundary  value  problems  can  be  expressed  as  variational 
integrals  to  be  minimized.  Accordingly,  one  possible  test  for  conver¬ 
gence  or  the  numerical  solution  to  the  boundary  value  problem  is 
whether  it  minimizes  the  associated  variational  integral. 


Problem 


The  primary  purpose  of  this  thesis  is  to  determine  the  usefulness 
of  the  variational  integral  as  a  test  for  terminating  l_he  iterative 
process. 

The  secondary  purpose  is  to  determine  the  usefulness  of  the 
variational  integral  as  a  means  of  determining  which  numerical  method 
gives  a  more  accurate  solution. 

Scope 

Tins  thesis  is  limited  to  the  study  of  stationary  (time-inde¬ 
pendent)  heat  equations  in  one  and  two  spacial  dimensions.  The  numeri¬ 
cal  approximations  to  the  boundary  value  problems  consist  of  the 
finite-difference  and  finite-element  techniques.  The  iterative  schane 
used  is  the  Gauss-Seidel  method  (method  of  successive  displacements) . 
Approach  and  Presentation 

Section  II  introduces  the  theory  necessary  to  understand  the 
variational  integral,  the  finite-difference  method,  and  the  finite- 
element  method.  The  Gauss-Seidel  iterative  technique  is  discussed. 

The  concept  of  norms  is  introduced  and  is  applied  in  examining  stopping 
criteria. 

Section  III  introduces  the  three  stationary  heat  equations  con¬ 
sidered  in  this  thesis.  The  variational  integral  formulations  of  the 
problems  are  given  and  the  numerical  approximation  techniques  are 
applied. 

Section  IV  presents  the  results  of  the  study.  First,  errors 
caused  by  approximating  the  variational  integral  are  examined.  Next, 
the  utility'  of  the  variational  integral  as  a  stopping  criterion  is 


discussed.  Finally,  the  usefulness  of  the  variational  integral  as  a 
means  of  determining  which  numerical  method  gives  a  more  accurate 
solution  is  analyzed. 

Section  V  draws  conclusions  and  gives  reccrrnenda ■cions  for  future 


studies. 


Theory 


Variations!  Integral 

I-Iany  problems  in  the  physical  sciences  can  be  solved  in  two  ways. 
First,  they  can  be  written  as  differential  equations  subject  to  given 
boundary  conditions.  Secondly,  they  can  be  written  as  variational 
integrals  to  be  extremized  (maximized  or  minimized) .  The  two  formula¬ 
tions  are  equivalent  because  the  function  that  will  extranize  the 
variational  integral  will  also  satisfy  the  differential  equation.  This 
equivalence  ia  demonstrated  by  the  fact  that  the  variational  integral 
is  extremized  only  when  an  equation,  called  the  Euler-Lagrange 
equation,  is  satisfied.  This  Euler-Lagrange  equation  is  precisely  the 
same  as  the  differential  equation  formulation  of  the  problem 
(Ref  8:67).  In  general,  only  elliptic  differential  equations  can  be 
written  in  variational  form  (Refs  8:117;  5:164)  but  seme  exceptions 
have  been  found  (Ref  7:252-256). 

In  variational  calculus,  one  looks  for  the  function  t(x,y)  with 
continuous  second  partial  derivatives,  that  extremizes  I(t)  on  region  R 
of  the  x-y  plane: 


where 


I(t)  =  //RF(x,y,t,tx,ty)dxdy 


t  =  t(x,y) 


(2.1) 


(2.2) 


(2.3) 

(2.4) 


To  find  t(x,y),  examine  the  set  of  functions 


t(x,y)  =  t (x,y)  +  en(x,y) 


(2.5) 


where  £  takes  on  sane  value  close  to  zero  and  n(x,y)  is  an  arbitrary 
function  with  continuous  second  partial  derivatives.  At  this  point,  an 
important  assumption  is  made.  It  is  assumed  that  t  is  given  on  the 
boundary  as  a  constant  or  as  a  function  of  x  and  y.  As  a  consequence, 
t  is  not  subject  to  variation  on  the  boundary  and  n(x,y)  must  equal 
zero  on  the  boundary. 

Substituting  Eq  (2.5)  into  Eq  (2.1)  gives 


J(e)  =  l(t+en)  =  //R] F(x,y, t+en,tx+enx,ty+eny)dxdy  (2.6) 

A  necessary  condition  for  I  to  ;iave  an  extremum  at  t  =  t  is  that 
J(e)  must  have  an  extremum  at  £  =  0.  So 


Id1™  e=G  =  0 


(2.7) 


Differentiating  Eq  (2.6)  with  respect  to  e  gives 


3  \  _  rr  /  3F  3t  ,  3F 

5?™  -  3l+  3T 

X 


dtx  _3F_ 
3e  3ty 


3t 

^rjdxdy 


(2.8) 


or 


9  T/  \ 


3t  X3tx  y3tyy 


(2.9) 


By  the  chain  rule  (Ref  1:92) 


therefore 


and 


3  /  3F  \  _  3  /  3F  \ 

9x(^j-  nax(atjc  J 

3  /  3F  \  3  /SF  \ 


Substituting  Bqs  (2.11)  and  (2.12)  into  Eq  (2.9)  gives 


By  applying  Green's  Theorem  (Ref  1:197) ,  the  second  integral  in  Eq 
(2.13)  can  be  transformed  into  a  line  integral 


5  $ 


Recall  that 


if  I  is  to  have  an  extremum  at  t  =  t.  Therefore 


(2.7) 


/;Rn  dxdy_  0 


(2.16) 


Since  n  is  an  arbitrary  function  of  x  and  y  then 


Eq  (2.17)  is  the  Euler-Lagrange  equation  and  it  must  be  satisfied  for  t 
to  extremize  I(t) .  In  the  derivation  that  lead  to  Eq  (2.17) ,  it  was 
assumed  that  the  function,  t,  was  specified  on  the  boundaries  (Dirich- 
let  boundary  conditions) .  Other  boundary  conditions  could  produce  a 
change  in  the  integral  that  is  to  be  extremized.  A  method  of  determining 
the  variational  integral  for  other  boundary  conditions  is  given  in 
Mikhlin  (Ref  12:116-121).  For  all  the  problems  considered  in  this 
thesis,  the  Euler-Lagrange  equation  given  as  Eq  (2.17)  will  give  us  a 
function,  F  that  extremizes  I(t)  as  given  in  Eq  (2.1). 

Finite-Element  liethod 

The  finite  element  method  involves  extremizing  the  variational 
formulation  of  the  problem.  The  solution  region  is  divided  into 
interconnected  subdomains  —  the  finite  elements.  The  elements  in 
one-dimensional  problems  are  line  segments.  The  most  versatile  ele¬ 
ments  in  higher  dimensions  are  triangles  for  two-dimensional  problems 
and  tetrahedrons  for  three-dimensional  problems  (Ref  8:86).  I  text, 
iterpolation  functions  are  chosen  to  represent  the  temperatures 


over  each  element.  These  interpolation  functions  are  often  poly¬ 
nomials.  Suppose  the  function  being  integrated  under  the  variational 
integral  contains  derivatives  up  to  the  (r+l)th  order.  Then,  the 
interpolation  functions  must  have  rth  order  continuous  derivatives  at 
the  element  interfaces  and  (r+l)th  order  continuous  derivatives  within 
each  element  (Ref  8:124).  When  these  criteria  are  met,  the  variational 
integral  for  the  entire  solution  region  can  be  written  as  the  sum  of 
the  variational  integrals  of  each  element  (Ref  8:78). 

An  interpolation  function  for  the  temperatures  on  an  element  is 
written  in  terms  of  the  nodal  values  associated  with  that  element.  The 
variational  integral  is  evaluated  over  the  element  and  the  integral  is 
differentiated  with  respect  to  the  nodes  and  set  equal  to  zero  in  order 
to  extremize  the  variational  integral  on  that  element.  The  equations 
for  each  element  are  assembled  to  give  the  equations  governing  the 
temperature  over  the  entire  region.  In  this  manner,  a  set  of  equations 
is  created  which  can  be  solved  simultaneously  to  find  the  temperatures 
at  the  element  nodes.  This  set  of  simultaneous  equations  can  be 
written  in  natrix  form 

At  =  b  (2.18) 

where 

A  =  a  coefficient  matrix 
t  =  vector  of  unknown  nodal  temperatures 
b  =  vector  of  known  constants 

Rather  than  go  into  the  details  of  the  process  here,  an  example  is 
given  in  Appendix  A. 
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Finite-Difference  Method 

A  differential  equation  can  be  converted  into  finite-difference 
tom  by  approximating  the  derivatives  in  the  equation  by  differences. 
The  elementary  definition  of  a  derivative  is  (Ref  3:54) 


dt  lim  t(x  +  Ax)  -  t(x) 


Ax 


Hie  derivative  is  approximated  by  a  divided  difference 


dt  ~  t(x  +  Ax)  -  t(x) 
dx  Ax 


(2.2C 


In  one-dimensional  problems ,  the  solution  region  is  broken  into 
line  segments  and  a  nodal  point  and  its  neighboring  nodal  points  are 
used  to  approximate  the  derivative  at  that  nodal  point.  In  this 
manner,  a  linear  difference  equation  is  created  for  each  nodal  point. 
The  equations  for  the  individual  nodal  points  can  be  assembled  to  give 
a  set  of  algebraic  equations  to  be  solved  simultaneously.  Linear 
combinations  of  Taylor  series  expansions  give  several  difference 
formulas  with  their  associated  truncation  error  (Ref  3:55-59): 

First  forward  difference 


dt  _  t (x  +  h)  -  t (x)  ,  nll_s 
dx - E - +  0(h) 


(2.21 


First  central  difference 


dt  _  t (x  +  h)  -  t (x  -  h)  .  2x 


m  i'  rl 


nv 


n  v 


First  backward  difference 


dt  _  t(x)  -  t(x  -  h)  ^ 

a? - e - +  0(h) 


(2.23) 


Second  forward  difference 


t(x  +  2h)  -  2t(x  +  h)  +  t (x) 


+  0(h) 


Second  central  difference 


(2.24) 


d~t  _  t(x  +  h)  -  2t(x)  +  t(x  -  h)  +  „  2^ 

dx2~  h2 

Second  backward  difference 

d2t  _  t(x)  -  2t(x  -  h)  +  t(x  -  2h)  , 

— - - - +  o(h) 

cbi  hz 


(2.25) 


(2.26) 


where  h  is  the  distance  between  nodal  points.  Truncation  error  is  a 
result  of  ignoring  higher  order  terms  in  the  Taylor  series.  The 
finite-difference  method  is  also  subject  to  scnething  called  round-off 
error.  Found-off  error  is  the  result  of  the  fact  that  a  ccnputer  can 
store  a  limited  number  of  significant  figures.  The  total  error  due  to 
both  truncation  and  round-off  error  will  be  referred  to  as  the  discreti¬ 
zation  error.  It  is  the  difference  between  the  analytical  solution  and 
the  finite-difference  solution  at  a  particular  nodal  point. 

For  two-dimensional  problems,  a  rectangular  mesh  is  imposed  over 
the  solution  region  in  order  to  find  approximate  temperatures  at  the 
nodal  point  of  the  mesh. 


Both  the  finite-element  and  finite-difference  methods  lead  to  a 

set  of  simultaneous  liner  algebraic  equations.  There  are  two  methods 

of  solving  for  the  unknowns  in  these  equations  —  direct  methods  and 

indirect  methods.  Direct  methods  lead  to  an  exact  answer,  excluding 

round-off  error,  after  a  finite  number  of  arithmetic  operations. 

Direct  methods  tend  to  be  time  consuming  since  each  element  in  the 

3 

coefficient  matrix  must  be  operated  upon.  In  general,  N  operations 
are  required  where  N  is  the  number  of  equations  (Refs  3:111;  5:209). 
Sane  examples  of  exact  methods  are  Gaussian  Elimination  and  the  Thotes 
method  (Ref  3:9-12;  44-48). 

Indirect  iterative  methods  require  an  infinite  number  of  opera¬ 
tions  to  solve  the  equations  exactly,  but  an  adequate  solution  can 


usually  be  achieved  in  a  finite  number  of  steps.  These  methods  usually 
require  less  than  operations,  especially  if  there. are  many  zeros  in 
the  coefficient  matrix  as  is  the  case  in  the  finite-element  and  finite- 
difference  equations.  Furthermore,  cumulative  round-off  error  does  not 
grow  as  in  the  direct  methods  (Ret  3; 111) . 

The  matrix  equation 


At  *  b 


(2.18) 


can  be  re-written 


t  =  Gt  +  c 


(2.27) 


An  iterative  solution  can  be  written  in  the  form 


t (n+1) 


(2.28) 


where 

t  ^  =  nth  iterated  vector 
G  =  iteration  matrix 
c  =  kncwn  constant  vector 

The  three  most  ccmonly  used  iteration  methods  are  the  Jacobi, 
Gauss-Seidel,  arid  successive  over-relaxation  methods.  They  are  derived 
by  factoring  the  coefficient  matrix  into  strictly  upper  triangular, 
diagonol  and  lcwer  triangular  matricies  denoted  by  U,  D,  and  L  respec¬ 
tively  (Ref  3:46-47): 


(U  +  D  +  L)  t  =  b  (2.29) 

Hie  Gauss-Seidel  toethod  (also  kncwn  as  the  method  of  successive 
displacements)  is  used  in  this  thesis.  It  expresses  the  (n+l)th 


12 


iterative  vector  in  terms  of  both  the  nth  iterative  vector  elements  and 
the  (n+1) th  iterative  vector  elements  as  soon  as  they  become  available: 


.  (n+1) 


=  -Ut' 


(2.30) 


or  written  in  a  form  similar  to  Eq  (2.27) 

t(n+1)  =  -(L  +  D)_1Utln)  +  (L  +  D)_1b 
Translated  into  a  computational  algorithm: 


(2.31) 


(2.32) 


where  a. _  is  the  element  in  the  rth  row  and  the  jth  column  of  the 
coefficient  matrix.  One  begins  the  process  by  choosing  seme  arbitrary 
initial  trial  vector,  t ^  . 

Norms 

In  the  following  sections,  we  will  be  concerned  with  measuring  the 
difference  between  exact  temperatures  derived  analytically  and  the 
temperatures  found  from  the  iterative  process.  The  error  vector  is 
defined  as 


(2.33) 


where 

e ^  =  the  error  vector  for  the  nth  iteration 
t  =  vector  of  exact  nodal  temperatures 

i  „  v 

t  r  =  vector  of  nodal  tempera tures  from  the  nth  iteration 
The  displacement  vector  (also  known  as  the  "residual  vector")  is 


defined  as 
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t'  -  vector  of  nodal  temperatures  from,  the  (n+l)th  iteration 
/  _  \ 

t'  1  =  vector  of  temperatures  frcm  the  nth  iteration 
It  is  often  more  convenient  to  measure  the  size  of  these  vectors 
as  scalar  quantities.  A  norm,  symbolized  by  |  |  •  |  | ,  will  allow  one  to 
do  just  that.  A  vector  n-norm  can  be  defined  try 

lltIL  -  £|t .|n)1/n  (2.35) 

i 

where 

t  =  vector 

t^  =  ith  element  of  vector  t 
The  norm  used  in  this  thesis  is 


lltllj  =  Ht.l  (2.36) 

i 

Stepping  Criteria 

Standard  Method.  One  of  the  fundamental  difficulties  with  itera¬ 
tive  methods  is  determining  when  to  step*  iterating.  Same  cannon 
criteria  are  (Fef  4:227) 


t<n>  .  t <  c 


(2.37) 


or 


fc(n)  _  t(n-l) 


.  (n) 


<  e 


(2.38) 


for  sane  prescribed  e  or 


n  >  N 


(2.39) 


for  sane  given  N. 

The  problem  with  the  above  criteria  is  that 


does  not  imply  that 


t  h-l) 


<  e 


-  t I)  <  e 


(2.37) 


(2.40) 


Figure  2.2  illustrates  this  fact  (Ref  9:5).  Thus,  for  problems  where 
convergence  is  slew,  tl  '  may  be  far  frem  the  exact  temperature  even 
though  it  would  have  eventually  been  reached.  However,  choosing  an 
iteration  limit  that  is  to  stringent  might  lead  to  a  considerable  waste 
of  computer  time. 

A  simple  technique  has  been  found  that  will  allow  one  to  use  the 
displacement  vector,  j  |d/n^  |  |  to  estimate  the  error,  [  je^  j  |  (Ref  14) . 
Recall  that  the  matrix  equation  can  be  written  as 


t  =  Gt  +  c 


(2.27) 


/nd  an  iterative  solution  can  be  written  as 


t(n+D  _  Gfc(n)  +  c 


(2.28) 


By  subtract  Eq  (2.28)  frcm  Eq  (2.27),  the  error  vector  can  be  written 


e(n+1)  =  Ge(n) 


(2.41) 


The  following  discussion  holds  for  all  real  symmetric  and  seme 
real  nonsyrmetric  matricies  (Ref  3:27).  Assume  that  the  initial  error 


/ 

V  +- 
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i 

i 

j 

i 
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iteration  limit 


Iteration  I, tracer 


Figure  2.2.  Iteration  Limit. 


vector  can  be  expanded  in  terns  of  eigenvectors  of  the  iteration 


matrix  G  (Fief  14:236) 


e‘°>  =  £  C.v.  (2. 

i*l  1-1 

where  Cf  are  scalars  and  m  is  the  dimension  of  the  square  coefficien 


matrix.  From  Eq  (2.41) 


iSjj 


I 


2 


;(n)  =  ^(0) 


where  Gn  is  G  raised  to  the  nth  power. 


(2. 


Blit  GVt  =  X  V,  by  the  definition  of  an  eigenvalue  where  X.  is  an 
eigenvalue  corresponding  to  V^.  Then 


m 


(n)  =  Z  C.X.’V. 


i=l 


li—i 


s<n)  -  hn{f& +i^ \n°&  *-+( 


(2. 


If  |  A^|  >  ! X2 1  i  1^3'  1  •••  1  I \n' 
then  for  large  sufficiently  large  n 


I'"’  *  hVl 


(2. 


so 


e 


(n+1)  _  ,  n+ln  .. 

~  Ai  LiXi 


e  (n+D 


In  order  for  the  error  vector  to  vanish 


(2. 


\\ |  <  1  (2. 

|  A^  |  is  called  the  spectral  radius.  Eq  (2.50)  shows  that  if  X1  is 
small,  convergence  occurs  rapidly.  In  general,  increases  for  a 
larger  number  of  unkncwns  (kef  3:116;  16:87).  Thus,  a  problem  takes 
nore  iterations  to  solve  as  one  increases  the  number  of  nodal  points 
used  in  the  finite-difference  and  finite-element  methods. 


Fran  Eq  (2.33) 


t  =  ttn)  +  e'n) 


and  fran  Eq  (2.50) 


_  .  (n+1)  (n) 

t  ~  t  +  X^e 


In) 


eliminating  e  between  Eq  (2.52)  and  Eq  (2.53)  gives 


t  = 


t(n+1)  -  aj t< n) 


l  -  x. 


(2.5: 


(2.53 


(2.54 


or 


(2.55 


From  Eq  (2.55) ,  it  can  be  seen  that,  if  A  is  close  to  one,  then  snail 
differences  between  successive  iterations  does  not  imply  that  the 
iterative  solution  is  close  to  the  exact  solution.  But,  given  the 
definition  ox  the  error  vector  and  the  displacement  vector,  the  error 
vector  can  be  approximated  by 


For  most  problems  A  will  not  be  kncwn  but  can  be  approximated  in 


or 


Therefore 


(Kef  14:241). 

For  sufficiently  large  n 

e(n)  = 

Xje'"-1' 

(2.58) 

(n+1)  (n)  * 

-  e 

(2.59) 

(n+1)  fc(n)  ~ 

(2.60) 

d(n)  * 

Xjd'1'-1' 

(2.61) 

i,  i  -  J. 

ld(n) 1 1 

(2.62) 

I'1'  ii 

d<n-wll 

|e(n)||  s  _ 

1  la(nl  1 1 

(2.63) 

IS.  »  ’ 

1 |a|n,l 1 \ 

l1 

!id(n_1)  !  ) 

One  new  has  a  simple  means  of  estimating  the  error  between  the 
true  solution  and  the  iterative  solution  at  the  nth  iteration  and, 
therefore,  a  criterion  for  stopping  the  iterative  process.  The 
variational  integral  will  be  judged  against  this  criterion.  Young 
provides  a  more  refined  version  of  this  error  approximation  stepping 
criterion  (Kef  17:226-227).  Warren  gives  a  more  detailed  analysis  of 
this  and  related  error  estimates  (Kef  16) . 

Variational  Integral.  Since  the  variational  integral  is 
extrema  zed  for  the  exact  solution  to  the  boundary  value  problem,  it  is 


.  V.  V  _ 
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thought  that  the  nodal  temperatures  from  the  iterative  nethod  can  be 
substituted  into  a  numerical  approximation  of  the  variational  integral. 
As  the  nodal  temperatures  approach  the  exact  values  at.  the  nodal 
points,  the  variational  integral  should  be  extremized.  Thus,  an  error 


criterion  might  be 


(n)  _  <  c 


(2.64) 


for  sane  prescribed  e.  This  criterion  might  not  have  the  difficulty 
that  the  standard  criteria  had  for  slcwly  converging  problems  as 
illustrated  in  Figure  2.3. 


Figure  2.3.  Hypothetical  Iteration  Limit  Comparison 


It  is  also  thought  that  the  variational  integral  can  be  used  as  a 
stopping  criteria  for  a  special  case  of  problems  where  the  temperatures 
iron  the  iterative  technique  are  either  mostly  above  or  mostly  below 
the  exact  values.  Suppose  the  initial  trial  solution  for  nodal  tempera¬ 
tures  is  above  the  exact  values  but  the  iterative  method  converges  to 
temperatures  below  the  exact  values;  then,  at  some  iteration,  the 
iterative  solution  will  pass  through  the  exact  solution  and  the  varia¬ 
tional  integral  will  be  extremized.  In  practice,  one  will  not  know  the 
exact  values  but  an  initial  trial  solution  can  be  chosen  that  is  known 
to  be  above  the  true  solution  and  then  an  initial  trial  solution  can  be 


chosen  that  is  known  to  be  be  lew  the  true  solution.  One  of  these  trial 
solutions  will  pass  through  the  exact  solution. 


Ill  Procedure 


In  this  section,  the  three  heat  equations  used  in  this  thesis  an< 
the  computer  codes  used  to  evaluate  the  problems  are  examined . 
One-Dimensional  Heat  Equation 

The  problem  to  be  examined  is  shown  in  Figure  3.1.  The  equation 
to  be  solved  is 

-  m2t  =  0  (3.1 

dx2 

where  m  is  a  constant.  The  boundary  conditions  are 


The  equation  is  easily  solved  (Ref  11:21): 

t  =  oosh[m(L  -  x) ) 
cos'  (mL) 


(3.4) 


Variational  Integral.  The  variational  integral  for  the 
one-dimensional  heat  equation  can  be  developed  by  comparing  the 
governing  differential  equation  (Eq  (3.1))  with  the  Euler-Lagrange 
equation  in  one  dimension: 

at  asht-/-0  ( 

The  result  of  the  comparison  is  that 

r-^2  +  H%)2 


Thus,  the  integral  to  be  extremized  can  be  written 


1  L 

Kt) 


(3.7) 


In  order  to  use  the  variational  integral  when  only  the  nodal 
temperatures  are  kncwn,  the  derivative  in  Eq  (3.7)  is  approximated  by  a 
forward  difference  and  the  integral  is  approximated  by  a  sum  (Ref 
10:235): 


1N 

I(t)  =  j  Z 
1=0 


\  xi+l  "  XJ 


.  2.  2 

+  m  t^ 


(xi+l  "  xi> 


(3.8) 


One  can  use  the  boundary  condition  given  in  Eq  (3.3)  and  the  forward 
difference  approximation  given  by  Eq  (2.21)  to  shew  that  for  i  =  N, 


For  the  purpose  of  validating  trie  computer  code,  the  variational 
integral  was  solved  analytically  in  Appendix  A: 

I (t)  =  tanh(mL)  (3.1( 


Finite-Difference  Forrrulation.  The  thin  rod  shown  in  Figure  3.1 
is  divided  with  a  number  of  equally  spaced  nodes  as  shown  in  Figure  3 


Figure  3.2.  Nodal  Divisions  for  Thin  Fod. 

The  differential  equation  is  approximated  using  Eq  (2.25): 


Uq  (3.11)  can  also  be  written  as 

"Vi  +  -  ti+i  =  0  {3-13) 

where 

D  =  2  +  (mhL) 2  (3.14) 

Using  the  boundary  condition  given  as  Eq  (3.3)  and  the  central 
difference  approximation  given  by  Eq  (2.22) , 

^J+l  =  *11-1 

thus  Eq  (3.13)  for  i  =  N  is  given  by 

"2Vl  +  ^  - 0 

and  using  the  boundary  condition  given  as  Eq  (3.2) 

Dt  -  t„  =  1  (3.17) 

1 

Therefore,  for  a  problem  with  four  unknown  nodal  temperatures ,  the 
simultaneous  equations  can  be  written  in  matrix  form  as  shewn  in  Figure 


D  -1  0  0 

-1  D  -1  0 

V 

fc2 

1 

0 

0  -1  D  -1 

0  0  -2  D 

fc3 

0 

L°_ 

Figure  3.3.  Matrix  Equation  for  Finite-Difference  Approximation,  Four 


(3.15) 


(3.16) 


Modes. 


In  order  to  validate  the  computer  code,  m  =  2  and  L  =  1  so 
that  the  temperatures  could  be  ccrpared  with  those  given  in  the 
literature  (Kef  11:243). 

Finite-Element  Formulation.  The  finite-element  method  is  more 
involved  than  the  finite-difference  method  and  the  derivation  of  the 
matrix  equation  is  given  in  Appendix  E.  The  matrix  formulation  of  the 
problem  turns  out  to  be  identical  to  the  finite-difference  formulation 
(Figure  3.3)  except  that 


„  _  12  +  4  (mhL) ' 

U  —  - *r 


6  -  (mhL)  ‘ 


(3.18 


Two-Dimensional  Poisson's  Ecuation 


The  problem  to  be  examined  is  a  square  plate  with  uniform  energy 
generation  across  the  plate.  The  equation  to  be  solved  is 


3x^  dyZ  K 


where  g  and  k  are  constants.  The  boundary  conditions  are 


|£<0,y)  =  0 


0<x,O> 


t(L,y)  = 


t  (x,L)  = 


The  analytic  solution  is  rather  involved  and  is  aeveloped  in  Appendix 
C.  For  the  validation  purposes,  g,  k,  and  L  were  set  equal  to  one  so 
that  the  numerical  solutions  could  be  compared  with  the  literature  (Ref 
11:260,384) . 

Variational  Integral.  The  variational  integral  is  developed  by 
comparing  the  differential  equation  with  the  Euler -Lagrange  equation  in 
two  dimensions  (Eq  2.17)).  The  integral  is  given  by 


The  integral  is  approximated  by  sums  over  both  the  x  and  y  directions. 
The  derivatives  are  approximated  by  backwards  differences  in  order  to 
take  advantage  of  the  derivative  boundary  conditions  in  the  same  way 
that  the  forward  difference  was  used  to  take  advantage  of  the  deriva¬ 
tive  boundary  condition  in  the  cne-dimensionai  problem  (pc.  23) : 


.  N  M 
I(t)  I 

i=0  j=0 


In  addition,  the  variational  integral  was  solved  analytically  for 
computer  code  validation.  The  derivation  is  very  similar  to  the 
derivation  for  the  one-dimensional  problem  given  in  Appendix  A. 

Finite-Difference  Formulation.  The  square  plate  is  overlayed  by  a 
system  of  square  meshes  each  of  dimension  h  by  h. 


27 


The  heat  equation  Eq  (3.19))  can  be  approximated  using  Eq  (2.25): 


i_(t 

h2lti+l,j 


-  2ti,j  +  4-1. j> 


¥fci,j+l  -  2ti,j  +  *1, j-l> 


+  |=0  (3.27) 


where 

x  =  ih  (3.28) 

y  =  jh  (3.29) 

The  subscripts  must  be  combined  before  the  function,  t,  can  be  repre¬ 
sented  by  a  vector.  This  is  easily  done  by  incrementing  tiie  subscript 
j  through  its  range  of  values  while,  for  each  j,  incrementing  i  through 
its  range  of  values.  A  matrix  equation  is  found  for  the  two-dimensional 
problem  in  a  manner  very  similar  to  that  of  the  one-dimensional  problem 
discussed  on  page  24. 

Finite-Element  Formulation.  Triangular  elements  are  used  in  the 
two-dimensional  problem.  If  one  uses  the  same  nodal  points  as  used  in 
the  finite-difference  scheme,  then  the  elements  can  be  oriented  in  two 
different  fashions.  The  orientation  in  Figure  3.4  will  be  referred  to 
as  "Case  One"  and  the  orientation  in  Figure  3.5  will  be  referred  to  as 
"Case  Two. "  The  two  arrangements  lead  to  different  matrix  equations 
and,  consequently,  different  solutions. 

IVo-Dimensional  Laplace's  Equation 

The  problem  to  be  considered  is  a  square  plate 

—■  +  =  0  (3.30) 

3x^  3y^ 
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v  j: 


Figure  3.5.  Finite-Element  Arrangement  (Case  Two). 


The  problem  is  a  rather  standard  boundary  value  problem.  Thererore, 
the  analytic  solution  will  not  be  derived  in  this  thesis.  The  solution 
is  given  by  (Eef  6:107): 


t(x,y)  =  Z 
n=l 


20o(l  -  oos(nTr))  sinh  (  j  sin  (  j 


nTrsinh(nn) 


(3.35) 


Variational  Integral.  Hie  variational  integral  is  developed  by 
ccrrparing  the  differential  equation  with  the  Luler-Lagrange  equation. 
The  integral  is  given  by 


I(t)  '  l"R  (h)  -  (H)  ^ 


(3.36) 


.  .  v \ 


The  integral  is  approximated  by  sums  over  both  the  x  and  y  directions. 
The  derivations  are  approximated  by  central  differences. 


Finite-Difference  Formulation .  The  square  plate  is  overlayed  by  a 
system  of  square  meshes  and  the  problem  is  solved  in  the  same  way  as 
the  Poisson  problem.  Due  to  time  constraints,  the  problem  was  not 
solved  using  finite-element  techniques. 

Computer  Codes 

All  codes  were  written  in  Fortran-77  for  a  DEC  Vax  11/780  ccnputer. 
The  ccnputer  programs  were  quite  similar  for  all  three  heat  equations. 

For  each  of  the  boundary  value  problems,  a  program  was  developed 
that  calculated  the  exact  temperatures  at  the  nodes  and  then  varied 
these  temperatures  by  a  given  percentage.  These  temperatures  were 
substituted  into  the  variational  integral  approximation  to  determine 
what  effect  errors  in  the  temperatures  had  on  the  value  of  the  varia¬ 
tional  integral. 

Another  code  was  written  for  each  problem  that  was  used  in  the 
cinalysis  of  the  variational  integral  as  a  stopping  criterion  and  as  an 
accuracy  criterion.  The  matrix  equations  from  the  finite-difference 
and  finite-element  methods  were  solved  iteratively  using  the  Gauss- 
Seidel  method.  The  exact  solution  at  the  nodal  points  was  calculated 
using  the  analytical  solution.  After  each  iteration,  the  difference 
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between  the  exact  and  iterative  temperatures  was  used  to  calculate  the 
percent  difference  between  the  exact  and  approximate  temperatures 


(3.38) 


The  difference  between  successive  iterations  was  used  to  calculate 
the  percent  change  in  temperature 


ia<n> 


t(n+l)  _  t(n) 


(n) 


(3.39) 


The  reason  for  these  new  relationships  is  that  when  one  actually 
loses  a  stopping  criteria,  the  percent  error  in  a  solution  is  usually 
more  useful  than  the  absolute  error.  This  also  allows  a  comparison 
between  the  standard  stopping  criterion  and  the  variational  integral 
criterion  by  eliminating  the  magnitudes  of  the  temperatures  and  varia¬ 
tional  integrals.  Note  that  Eqs  (3.38)  and  (3.39)  are  very'  similar  to 
the  error  vector  and  displacement  vector  previously  defined.  In  the 
following  sections  the  terms  "error  norm"  and  "displacement  norm"  will 
refer  to  Eqs  (3.38)  and  (3.39).  The  relationship 


ild(n)ll 

| !d(n-1) | | 


(3.40) 


is  assumed  to  hold  although  it  cannot  be  rigorously  developed  as  was  Dq 


The  variational  integral  was  computed  after  each  iteration  by 
substituting  the  temperatures  from  the  iterative  process  into  the 
variational  integral  approximation.  To  compare  the  variational  inte¬ 
gral  as  a  stopping  criterion  to  the  stopping  criterion  defined  in  Eq 
(3.40) ,  the  percent  change  in  the  variational  integral  is  defined  as 


IV  Numerical  Results 


Variational  Integral  liinimum 

Recall  that  to  calculate  the  variational  integral  based  on  the 
nodal  temperatures,  the  integral  had  to  be  approximated  using  differ¬ 
ences  for  the  derivatives  and  sums  for  the  integrals.  These  approxima¬ 
tions  lead  to  the  concern  that  the  variational  integral  approximation 
might  not  be  extremized  at  the  exact  nodal  temperatures. 

This  was  found  to  be  the  case  for  all  three  problems  considered. 
The  exact  temperature  for  each  node  was  varied  by  a  given  percentage  of 
the  exact  temperature  and  the  variational  integral  was  calculated. 
Figures  4.1,  4.2,  and  4.3  illustrate  that  the  integral  for  the  one- 
dimensional  heat  equation  had  a  minimum  for  sane  given  values  of 
temperature  and  that  this  minimum  occurred  closer  to  the  exact  tempera¬ 
tures  as  the  number  of  nodes  increased.  Figures  4.4  and  4.5  illustrate 
that,  for  the  two-dimensional  problems,  the  variational  integrals  are 
minimized  and  these  minima  occur  at  temperatures  other  than  the  exact 
values.  Table  I  summarizes  the  results  for  all  three  problems.  For 
the  purpose  of  this  thesis,  one  variational  integral  approximation  will 
be  said  to  be  more  accurate  than  another  if  the  niinimum  in  the  varia¬ 
tional  integral  occurs  for  temperatures  closer  to  the  exact  values. 
Thus,  Table  I  shews  that  the  variational  integral  approximation  becomes 
more  accurate  as  the  number  of  nodes  increase.  In  practice,  one  will 
not  knew  the  exaert  tereperatures  and  will  have  no  way  of  determining  the 
accuracy  of  an  integral  approximation  as  was  done  in  Table  I. 


The  exact  value  of  the  variational  integral  is  calculated  using 
the  analytical  solution  as  shewn  in  Appendix  A.  It  can  be  seen,  in 


Table  I,  that  the  value  of  the  variational  integral  approximation  does 
not  approach  the  value  of  the  exact  variational  integral  in  the  Poisson 
problem  even  though  the  accuracy  of  the  integral  approximation  improves 
for  an  increasing  number  of  nodes  (using  the  definition  of  accuracy 
given  above) .  Thus,  relating  the  value  of  the  integral  approximation 
to  the  exact  value  is  of  little  use.  Furthermore,  one  will  not  be  able 
to  calculate  the  exact  variational  integral  in  "real-world"  problems 
since  the  analytic  solution  will  not  be  knewn. 


Table  I 

Variational  Integral  Approximation 


Problem 

Exact  . 
Integral1 2 3 

Nodes 

Integral  2 

Approximation 

Percent  , 
Deviation1 

4 

1.260 

-0.80 

One- Dimen. 

0.964 

8 

1.105 

-0.21 

16 

1.033 

-0.05 

32 

0.993 

-0.01 

IW- Dimen. 

4x4 

-0.1070 

+13 

Poisson 

-0.1406 

6x6 

-0.0947 

+9 

8x8 

-0.0884 

+7 

IVo-Dimen. 

A 

4x4 

4.35x10, 

-53 

Laplace 

1.04x10^ 

6x6 

5.92x10, 

-45 

8x8 

7.21xl0J 

-39 

1.  Value  of  variational  integral  using  analytic  solution. 

2.  Value  of  variational  integral  approximation  at  its  minimum. 

3.  Percent  deviation  from  exact  nodal  temperatures  where  the 
minimum  in  the  variational  integral  occurs. 


Criteria 


Recall  from  Section  II,  that  for  a  problem  with  slew  convergence, 
the  difference  between  temperatures  for  successive  iterations  could  be 
small  even  though  the  difference  between  the  iterated  value  and  the 
exact  value  was  still  quite  large.  It  was  thought  that  the  variational 
integral  might  give  an  alternative  evaluation  for  termination  since  the 
integral  is  extranized  as  the  iterated  nodal  temperatures  approach  the 
exact  nodal  temperatures .  Figure  4.6  shews  that  as  the  sum  of  the 
iterated  temperatures  approaches  the  sum  of  the  exact  temperatures,  the 
variational  integral  converges  to  a  constant  value. 

Figure  4.7  compares  a  ccrrmon  stopping  criterion,  I  !d^l  I  with  the 
variational  integral  stopping  criterion,  Ul  n; !.  It  can  be  seen  from 
Figure  4.7  that,  past  a  certain  iteration,  |  |d*n^|  |  and  |ul^  | 


decrease  at  the  same  rate. 

One  can  easily  explain  the  linear  decrease  in  ! !d* 


(keeping  in 


mind  that  the  y-axis  in  Figure  4.7  is  a  logarithmic  scale) .  The  slope 
of  !  Id^  |  |  is  given  by 


,  (n-1) 


n  -  (n-1) 


(4.1) 


Recall  that  for  large  n 


ld(n)ll 

I j  (n-1) 


(4.2) 


log | A x [  =  log||d(n)j|  -  log] |d 


(n-1) 


(4.3) 


Iteration  '  Imt-er 


By  comparison  o£  Lq.  (4.1)  and  (4.3),  it  can  be  seen  that 


ni  =  log  | 


(4.4) 


More  importantly,  Figure  4.7  shews  that  for  sufficiently  large  n 


UT(n)  I 

x  |  ~  — — L_ 

l1  | , j (n-l) 


(4.5) 


Eq.  (4.5)  was  found  to  hold  lor  all  three  heat  problems  and  all  nodal 
densities  considered.  Thus,  the  variational  integral  does  not  provide 
a  better  stopping  criterion  because  |AI^  !  decreases  at  the  same  rate 
as  I  |d(n) | | . 

Furthermore ,  the  standard  method  can  provide  an  estimate  of  the 
error  while  the  variational  integral  cannot.  Figure  4.8  shews  the 
error  norm  and  the  error  norm  estimate.  The  error  estimate  approxi¬ 
mates  the  error  once  the  spectral  radius  approximation  has  converged  to 
its  final  value.  Figures  4.9  through  4.11  snow  that  the  spectral 
radius  approximation  rapidly  reaches  a  given  value  and  that  this  value 
increases  as  the  number  of  nodes  increase,  returning  to  Figure  4.8,  it 
can  be  seen  that  the  estimate  of  the  error  norm  continues  to  decrease 
even  though  the  error  norm  has  stopped  decreasing.  This  is  due  to  the 
fact  tliat  the  error  norm  estimate  depends  only  on  temperatures  from 
successive  iterations.  The  error  norm  decreases  until,  at  some  itera¬ 
tion,  the  discretization  error  prevents  any  further  decrease. 

In  Section  II,  it  was  predicted  that  the  variational  integral 
could  be  useful  when  the  iterative  temperatures  passed  through  the 
exact  values.  Figure  4.12  shows  the  sum  of  the  nodal  temperatures  fretn 


Iteration  i\'ur.l>er 


4.10.  Gnectral  Radius  Approximation,  1-1)  Problem,  Ripld:  1  lodes 


uoiriuro: 


the  iterative  method  passing  through  the  sum  of  the  exact  nodal  temper¬ 
atures  after  thirty-seven  iterations.  Thus,  one  would  expect  the 
variational  integral  to  be  extnemized  at  the  thirty-seventh  iteration. 
Figure  4.13  shews  that  the  variational  integral  was  actually  minimized 
after  nineteen  iterations.  Tins  is  understandable  in  light  of  the 
previous  discussion  in  which  it  was  found  that  the  variational  integral 
approximation  is  minimized  for  temperatures  other  them  the  exact 
values.  The  results  from  Figures  4.12  and  4.13  can  be  represented  more 
concisely  using  |AI^  |  and  |  |e^  [  J.  Figure  4.14  shows  that  the 
variational  integral  stopping  criterion,  |AI^|  has  a  "cusp"  or  local 
minimum  where  the  variational  integral  is  a  minimum  (nineteenth  itera- 

~  /n  \ 

tion)  and  the  error  norm,  |  |e'  '  |  |  has  a  cusp  where  the  iterative 
solution  passes  through  the  exact  solution  ( thirty-seventh  iteration) . 
In  Figure  4.14,  the  initial  trial  solution  was  greater  than  the  exact 
solution.  The  finite-difference  equations  converged  to  temperatures 
belcw  the  exact  values.  Figure  4.15  shews  that  when  the  initial  trial 
solution  was  belcw  the  exact  solution,  the  iterative  temperatures  did 
not  pass  through  the  exact  values  and  no  cusps  were  observed  for  either 
|£l(n)  |  or  |  |e^  |  | .  Table  II  ccnpares  the  iteration  number  where  the 
integral  approximation  is  minimized  to  the  iteration  number  where  the 
error  norm  is  minimized.  Appendix  D  shews  that  the  results  of  the 
stepping  criteria  study  hold  for  all  three  problems  and  all  nodal 


densities  considered. 


6 

15 

26 

One-Dimen. 

8 

30 

53 

10 

51 

90 

16 

150 

280 

TWo-Dirren. 

4x4 

20 

38 

Poisson 

6x6 

43 

96 

8x9 

74 

185 

Two- Dimen. 

4x4 

7 

32 

Laplace 

6x6 

18 

66 

8x8 

35 

111 

1.  Using  the  finite-difference  method 

2.  Iteration  number  where  the  minirnzn 

occurs 

Accuracy  Criterion 

Since  the  variational  integral  is  minimized  for  the  exact  solu¬ 
tion,  it  was  thought  that  the  variational  integral  approximation 
could  be  used  to  determine  which  numerical  method  would  give  the  best 
approximation  to  the  true  solution.  Two  methods  were  considered — the 
finite  difference  method  and  the  finite -element  method.  Due  to  time 


constraints,  only  the  one-dimensional  problem  and  the  Poisson  problem 
were  considered.  Table  III  shews  the  results  for  the  one-dimensional 
problem. 


liable  III 

Integral  and  Error  Norm  for  1-D  Problem 


Nodes 


Method 


Integral 


3 

FD 

1.3720 

0.0759 

FE 

1.3695 

0.0812 

4 

FD 

1.2606 

0.0549 

FE 

1.2594 

0.0569 

5 

FD 

1.1967 

0.0429 

FE 

1.1960 

0.0436 

6 

FD 

1.1553 

0.0351 

FE 

1.1550 

0.0352 

7 

FD 

1.1264 

0.0298 

FE 

1.1262 

0.0295 

8 

FD 

1.1051 

0.0259 

FE 

1.1049 

0.0253 

1.  FD: 

Finite-Difference  Method 

FE: 

Finite-Element  Method 

Note  that  tor  nodal  densities  of  less  than  six,  the  temperatures 
from  the  finite-element  method  minimize  the  variational  integral  even 
though  the  finite-difference  method  minimizes  the  error  norm.  Nodal 
densities  greater  than  eight  were  considered  but  they  gave  the  same 
results  as  densities  six  through  eight  and,  consequently,  were  not 
included  in  the  table.  The  discrepancy  between  the  integral  and  error 
norm  for  nodal  densities  less  than  six  can  be  explained.  Table  TV 
shews  that  the  temperatures  for  the  finite-element  method  are  less  than 
the  exact  temperatures  while  the  temperatures  from  the  finite-dif¬ 
ference  method  are  greater  than  the  exact  temperatures. 


'fable  iv 

Nodal  Temperatures  for  1.-D  Problem  (Four  Modes) 


Node 

Exact 

Finite-Diff . 

Finite-El em. 

1 

0.6252 

0.6289 

0.6215 

2 

0.4102 

0.4150 

0.4051 

3 

0.2998 

0.3049 

0.2943 

4 

0.2658 

0.2710 

0.2604 

1 

Figure  4.16  is  based  on  Figures  4.1  through  4.5  and  shews  hew  the 
variational  integral  could  fail  to  predict  which  method  is  most  accu¬ 
rate.  For  nodal  densities  greater  than  or  equal  to  six,  the  finite- 
element  method  gives  more  accurate  temperatures  and,  therefore,  the 
variational  integral  approximation  correctly  predicts  the  more  accurate 
method. 

The  results  for  the  Poisson  problem  are  summarized  in  Table  V.  It 
was  found  that  one  finite-element  case  gave  the  most  accurate  tempera¬ 
tures  while  the  other  finite-element  case  gave  the  least  accurate 
answers.  For  all  nodal  densities  considered,  the  variational  integral 
correctly  predicted  which  method  gave  the  most  accurate  solution.  The 
difficulty  that  occurred  in  the  one-dimensional  problem  did  not  occur 
for  the  Poisson  problem  because  the  majority  of  temperatures  for  all 
three  methods  gave  temperatures  belcw  the  exact  values. 


Schematic  Illustration  of  Accuracy  Criterion  Failure 


T&ble  V 


Integral  and  Error  Norn  for  Poisson  Problem 


|  Nodes 


Method 


Integral 


2x2 


FE-2  -0.1632 

FD  -0.1680 

FE-1  -0.1710 


0.4147 

0.1891 

0.0844 


4x4 


8x8 


FE-2  -0.1806  0.3966 

FD  -0.1815  0.2143 

FE-1  -0.1821  0.1024 


FE-2  -0.1940  0.3961 

FD  -0.1942  0.2366 

FE-1  -0.1944  0.1223 


1.  FD:  Finite-Difference  Method 

FE-1:  Finite-Element  Method  (Case  One) 
FE-2:  Finite-Element  Method  (Case  IVvo) 


Conclusions  and  Reccoiendations 


Conclusions 

The  variational  integral  approximation  for  each  of  the  three 
boundary'  value  problems  was  minimized  for  temperatures  other  tlian  the 
exact  values.  As  the  number  of  nodes  was  increased,  the  approximations 
were  minimized  at  temperatures  closer  to  the  exact  values. 

As  a  stopping  criterion,  the  variational  integral  was  found  to  be 
less  effective  than  the  error  estimate  criterion.  For  a  sufficiently 
large  number  of  iterations,  the  variational  integral  stopping  criterion 
decreased  at  the  same  rate  as  the  error  norm  estimate  criterion.  But 
the  error  norm  estimate  provides  an  estimate  of  the  actual  error  norm 
for  a  limited  number  of  iterations  while  no  such  error  estimate  is 
available  iron  the  variational  integral  stopping  criterion.  Further¬ 
more,  the  error  norm  estimate  is  easier  to  calculate. 

For  the  special  case  where  the  temperatures  frcn  the  iterative 
methods  passed  through  the  exact  temperatures ,  it  was  found  that  the 
error  norm  was  minimized  at  a  different  iteration  than  was  the  varia¬ 
tional  integral  approximation .  This  was  expected  since  the  variational 
integral  approximations  were  not  minimized  for  the  exact  temperatures . 
Thus,  the  utility  of  the  variational  integral  as  a  stopping  criterion 
for  this  special  case  depends  on  the  accuracy  of  the  variational 
integral  approximation. 

The  variational  integral  approximation  did  not  correctly  predict 
which  method  gave  the  most  accurate  solution  in  .one  of  the  cases 
studied.  It  is  believed  that  this  failure  is  due  to  the  variational 


integral  approximation  being  minimized  at  temperatures  other  tliar,  the 
exact  values. 

Kecarnendations 

Many  of  the  difficulties  of  using  the  variational  integral  as  a 
stopping  criterion  or  an  accuracy  criterion  are  believed  to  be  due  bo 
the  fact  that  the  variational  integral  approximation  is  minimized  at 
temperatures  other  than  the  exact  values.  Thus,  attention  should  be 
focused  on  reducing  or,  at  least,  predicting  this  error.  Better 
methods  of  approximating  the  derivatives  and  integrals  under  the 
variational  integral  my  exist. 

For  the  accuracy  criterion,  consideration  should  be  given  to  using 
linear  interpolation  functions  to  transform  the  nodal  temperatures  into 
a  piecewise  continuous  function.  These  functions  can  be  substituted 
into  the  variational  integral  and  the  integral  can  be  solved  analyti¬ 
cally.  Other  techniques  that  solve  the  boundary  value  problem  over  the 
entire  solution  region  such  as  the  Fayleigh-Fdtz  and  Galerkin  methods 
can  also  be  substituted  into  the  variational  integral  to  determine 
which  method  gives  the  most  accurate  solution. 
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Appendix  A 

Exact  Variational  Integral 


The  variational  integral  for  the  one-dimensional  heat  equation  is 


I(t) 


=  i  A 

2  'o 


dtr  +  r2t2 


d xl 


dx 


(A.l) 


Integrating  the  first  term  by  parts  allcws  Eq  (A.l)  to  be  written 


I(t)  =  ~ 


dx 


L 

x=0 


1  /i 

2  x=0 


'dx 


(§) 


-m2ta 


(A.2) 


Since  the  heat  equation  is  given  by 


0 


(A.  3) 


Eq  (A.2)  reduces  to 


I(t)  =  f 


dt 

bd£ 


L 

x=0 


(A.  4) 


Using  the  derivative  boundary  condition 


=  0 


(A.  51 


the  integral  can  be  written  as 


Consequently 


(e)  =  xjfci  ~  Vj  +  fcj  I  fci 

xj  '  xi  xj  ~  xi 


dt  _  fcj  ~  fci 
dx  xj  "  xi 

The  integral  for  a  given  element  can  be  written 


(B.  3) 


(B.4) 


!<e>  = 

2  Xi 


+  J  /Xj  m2  ( Xj  -1-— itj-  +  -i - -x  ]  dx 

2  Xi  Xi  "  Xi  Xi  -  Xi  i 


(B.5) 


■(e)  : 


After  I  is  integrated ,  it  will  be  a  function  of  t^  and  t^  only.  It 
(e) 

minimize  I  ,  it  must  be  differentiated  with  respect  to  both  t^  and 
t y  One  can  arbitrarily  do  the  differentiation  first  and  then  the 
integration: 


3I(e)  *  1*1-^  ,/l  m2  2, 

3ti  (xj-x,)2  Xi  Uj  -  xt)2  '  3  1  l 


+  (xitj  -  2Xjti  +  xitj)x  +  (tA  -  tj)x  ] 


(B.6) 


The  integration  is  carried  out  to  give 


3l(e)  =  fci  ~  fcj 
9ti  xj  -  xi 


-  x^' 


+  iv  (xj 


3XjX^  +  3XjX2  -  x2l  (B. 7) 


which  can  be  simplified: 


31 


(e) 


at. 


x .  -  x. 
1  1 


(xj  -  xi> 


(2ti  +  tj) 


(B.  8) 


Similarly 


(tf  +  2tj> 


(B.9) 


Assembling  the  integral  of  the  elements,  we  find 

E  (e) 

I(tl'  t2'  -V  -  V  (B-10> 

e=l 

To  find  the  minimum  of  I,  differentiate  with  respect  to  each  m  and  set 
each  derivative  equal  to  zero.  Since  there  are  M  nodes,  there  will  be 
M  equations.  Each  equations  will  of  the  form 


31 


(1) 


31 


(2) 


31 


(a) 


31 


(b) 


31 


(E) 


3t_ 


m 


at 


m 


at 


m 


at 


m 


3t 


=  0 


(b.ii: 


m 


The  only  two  elements  that  will  contain  t  in  their  integrals  are  the 
elements  immediately  adjacent  to  node  m.  If  the  elements  denoted  (a) 
and  (b)  are  the  elements  on  either  side  of  node  m,  then  Eq  (B.ll) 
reduces  to 


31 


(a) 


31 


(b) 


9tm 


at 


=  o 


(B.  12] 


m 


For  element  (a) ,  use  Eq  (B.9)  and  let  i  =  m  -  1  and  j  =  m 


31 


(a) 


at 


m 


t  -  t  2 

m  m-1  ,  m  Ax  x 

- x - +  — z —  (t.  .  +  2t  ) 

Ax  6  m-1  m 


(B. 13] 


67 


For  element  (b) ,  use  Eq  (B.8)  and  let  i  =  n,  and  j  =  m  +  1 


31^.  =  +  bQx  2t 

3t  Ax  6  m  m+1 

m 


where 


Ax  =  x .  -  x. 
3  i 


Substituting  Eqs  (B.14)  and  (B.13)  into  (B.12)  gives 


where 


-t  ,+Dt  -  t  ^,=0 
m  -  1  m  m+1 


^  _  12  +  4m2  (Ax) 2 

U  2  2 
6  -  in  (Ax) L 


Eq  (B.16)  holds  for  all  interior  nodes.  For  the  first  node,  the 
boundary  condition  can  be  used: 


tQ  =  1 


For  node  M,  the  temperature  appears  only  in  element  (E) .  Setting  i  = 
M  -  1  and  j  =  M  and  using  Lq  (B.9) 


fEI  _  tK  ~  t::-i  n2ix  [t  .  , 

r - s — +  sr  +  2V 


-2tH  -  1  +  ‘h-!  “  0 

;mere  D  is  as  in  Eq  (B.17) .  One  new  has  the  M  linear  algebraic 
equations  to  solve  simultaneously. 


Appendix  C 


Derivation  of  an  Analytical  Solution  for  Poisson's  Equation 
The  equation  to  be  solved  is 


(C.l 

t(x,L)  =  0 

(C.2) 

t(L,y)  =  0 

(C.3) 

|i(0,y)  =  o 

(C.4 

-  ° 

(C.5 

Assume  that  the  solution  can  be  written  as  t(x,y)  =  v(x,y)  + 
w(x,y)  where  v  is  the  particular  solution  of  Poisson's  Equation  and  w 
is  the  solution  of  the  associated  homogeneous  equation  (Ref  16:239): 


(x,  0) 


W-w(x,0) 

dy 


Assume  v  is  of  the  form 

v(x,y)  =  A  +  Bx  +  Cy  +  Dx^  i-  Exy  + 
Substituting  Eq  (C.12)  into  Eq  (C.6)  gives 


Let 


and 


2D  +  2F  =  -  J 


F  = 


g 

k 


D  =  0 


The  other  coefficients  are  arbitrary  so  let 


v(x,y) 


=  2^_ 


2k  2k 


9  2 


so  that  v  reduces  to  zero  on  the  sides  y  =  0  and  y 
New  solve  for  w  using  separation  of  variables 

w(x,y)  =  X (x) Y (y) 


This  gives 


x'lx)  -  XX (x)  =  0 

/  T2  2  \  . 
X(L)  =  t-S-U, 


I 

Fi¬ 


fe.  19) 


Y~(y)  +  AY  (y)  =  0 


Y(L)  =  0 


Y'(L)  =  D 


Try 


Y  =  Acos(/X~  yj+  Bsin(>/X  yj 


Apply  boundary  conditions: 


Y"(0)  =  0 


therefore 


B  =  0 


and 


Y(L)  =  0 


consequently 


A 


( (2n  +  1) 


-)2 

2L 


n  =  0,  1,  2 


and 


tram  Eq  (C.15) 


wnere 


[  (2n  +  1)^-] 
2L 


(C.32) 


Applying  boundary  condition 


Therefore 


X'(0)  =  0 


X(x)  =  C[en:  +  e-rx] 


(C.33) 


(C.34) 


7TV 

X(x)  =  cosh[(2n  + 


(C.35) 


By  superposition  (Ref  15:7) 


w(x,y) 


Z  C  cos  (2n  +  l)^r  oosh  ( 

Oil  /.Li 

L  J  l 


(2n  +  1) : 


(C. 36) 


Z^pplying  the  last  boundary  condition  (Eq  (C.19)): 


t2  2 

gL  _  gy  _ 

2k  2k 


Z  Cnoos  (2n  +  °°sh  (2n  +  Dj 

n=0  L  L 


(C.  37) 


Multiplying  both  sides  of  Lq  (C.37)  by  the  term  cos  [(2m  +  l)^-] 
and  then  integrating  with  respect  to  y; 


fo  2F  003  (2m  +  1)§;  if  005  _(2m 
=  /0Ln=0Cnaos  (2m+1)l  005  [(2n+1)S 


ro  fr  °°s  fi2™ + 


cosh  (2n+l)2  dy 


(C.  38) 


For  m  ?  n,  the  right-hand  side  equals  zero.  Therefore 


r  =  gL^ _ (-l)n+1  32 _ 

z  (2n  +  1)  ii  oosh[  (2n  +  1)£] 

and 

T  2  „  (-l)n+1oos[(2n+l)^]cosh[(2n+l)^] 

w(x,y)  =  z  - 3 - ^ - -  2L 

n=0  (2n+l) ^oosh [ (2n+l)  j] 

Since 


t(x,y)  =  v(x,y)  +  w(x,y) 


theii 

t2  2 
t(x,y)  =  2L_  _  gy_ 
'  ,y'  2k  2R 


2 


oo 


(-1)  n+1oos  [  (2n+l) cosh  [  (2n+l)  ^1 


Appendix  D 
Stopping  Criteria 


The  following  figures  show  the  error  norm,  the  error  norm  esti¬ 
mate,  and  the  variational  integral  stopping  criterion  for  all  three 

heat  equations  solved  by  the  finite-difference  method.  The  figures 

* 

shew  that  for  a  sufficiently  large  number  of  iterations,  the  error  norm 
approximation  and  the  variational  integral  stopping  criterion  decrease 
at  the  same  rate. 

The  figures  also  show  that  when  the  iterative  solution  passes 
through  the  exact  solution,  the  variational  integral  achieves  a  local 
minimum  at  a  different  iteration  than  does  the  error  norm.  In  the  case 
of  Laplace's  problem,  the  error  norm  achieves  a  local  minimum  but  it  is 
not  pronounced  enough  to  be  visible  on  the  graphs. 
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